
#ifndef GRID_H
#define GRID_H

#include "stdafx.hpp"
#include "array2d.hpp"
#include <string>
#include <ostream>
#include <fstream>
#include <cstdlib>
#include <limits>
#include <cmath>
#ifdef CAVITY_USE_BOOST_RANDOM
    #include <boost/random.hpp>
#endif
#ifdef CAVITY_MODES_VECTOR_USE_PTR_VECTOR
    #include <boost/ptr_container/ptr_vector.hpp>
#endif

//Easy to change the types!
#ifndef CAVITY_GRID_DATA_TYPE
	#define CAVITY_GRID_DATA_TYPE long double
#endif
#ifndef CAVITY_GRID_SIZE_TYPE
	#define CAVITY_GRID_SIZE_TYPE uword
#endif

/*template<typename T, typename uint>
class basic_func2d
{
    public:
	virtual void operator()(uint x, uint y, array2d<T, uint>& g, T Kd = 1.00) = 0;
};*/

//typedef basic_func2d<CAVITY_GRID_DATA_TYPE, CAVITY_GRID_SIZE_TYPE> func2d;

//NOTE: basic_grid not yet used. I just wrote it as part of broad changes.
template< typename container_t,
		  typename data_t = typename container_t::data_type,
		  typename size_t = typename container_t::size_type >
class basic_grid
{
	const size_t xstride;
	const size_t ystride;
	const size_t h;
	const size_t w;
	container_t & container;
public:
	typedef container_t container_type;
	typedef data_t data_type;
	typedef size_t size_type;

	basic_grid(container_t & container_, size_t xstride_, size_t ystride_)
		: xstride(xstride_), ystride(ystride_), h(container_.height()/ystride_),
		  w(container_.width()/xstride_), container(container_)
	{ /*empty*/ }
	basic_grid(container_t & container_, size_t stride_)
		: xstride(stride_), ystride(stride_), h(container_.height()/stride_),
		  w(container_.width()/stride_), container(container_)
	{ /*empty*/ }
	data_t & operator()(size_t x, size_t y)
	{
		return container(x*xstride, y*ystride);
	}
	const data_t & operator()(size_t x, size_t y) const
	{
		return container(x*xstride, y*ystride);
	}
	const size_t height() const { return h; }
	const size_t width() const { return w; }
};

class grid : public array2d<CAVITY_GRID_DATA_TYPE, CAVITY_GRID_SIZE_TYPE>
{
  public:
	typedef array2d<grid::data_type,grid::size_type> base_type;
	typedef std::pair<grid::size_type,grid::size_type> point2d;
	grid(size_type w_in, size_type h_in, data_type fill = 0.0);
	static const unsigned int precision = std::numeric_limits<grid::data_type>::digits10;

	struct span;
	typedef std::vector<span> region;
	typedef std::pair< point2d, std::vector<point2d> > normal_pair;
    typedef std::vector<normal_pair> boundary;

	template<typename func2d>
	bool iterate_region(const region & reg, func2d & f);

	const static bool binary = 1;
	const static bool text = 0;
};

template<typename T1, typename T2>
grid::point2d make_point(const T1 & a, const T2 & b)
{
    return grid::point2d(a,b);
}

void load(grid & g, const std::string & path, bool textorbinary = grid::text);

void write_file(grid & g, const std::string & path, bool textorbinary, grid::region & reg);

struct grid::span
{
    grid::size_type x, y, dy;
    grid::data_type Kd;
    span(grid::size_type x_in, grid::size_type y_in, grid::size_type dy_in, grid::data_type Kd_in = 1.0)
            : x(x_in), y(y_in), dy(dy_in), Kd(Kd_in) {}
};

inline unsigned long thing_to_count(grid::boundary::const_reference ref) { return ref.second.size(); }

inline unsigned long thing_to_count(grid::region::const_reference ref) { return ref.dy; }

template<class vec_T>
grid::size_type num_points(const vec_T & v, unsigned int until = std::numeric_limits<unsigned int>::max())
{
    unsigned long int points = 0;
	for(unsigned int i = 0; i < v.size() && i < until; ++i)
		points += thing_to_count(v[i]);
	return points;
}

template<class vec_T>
grid::size_type first_more_than(const vec_T & v, unsigned long int quota, grid::size_type begin_index = 0)
{
	unsigned long int count = 0;
	for(unsigned int i = begin_index; i < v.size(); ++i)
	{
		count += thing_to_count(v[i]);
		if(count > quota)
			return i;
	}
	return v.size(); //does returning this here make sense?
}

template<typename func2d>
inline bool grid::iterate_region(const region & reg, func2d & f)
{
	for(unsigned int i = 0; i < reg.size(); ++i)	//over all spans
	{
		for(unsigned int j = 0; j < reg[i].dy; ++j)	//for each span
		{
			f(reg[i].x, reg[i].y + j, *this, reg[i].Kd);
			//	x,			y,		  grid,	dielectric constant
		}
	}
	return true;
}

grid & operator *= (grid & a, const grid::data_type & c);

#ifdef CAVITY_MODES_VECTOR_USE_PTR_VECTOR
   typedef boost::ptr_vector<grid> vec_of_grids;
   //So vector operations won't copy grid data
#else
    typedef std::vector<grid> vec_of_grids;
    //grids will be copied like any other variable in a vector
#endif

// THE FUNCTION OBJECTS
/////////////////////////////////////////////////
// These are handy little things that can be passed to
// a grid's iterate functions.
////////////////////////////////////////////////

namespace functions
{

class relax //: public func2d
{
public:
	relax(grid & lap_ref, grid::data_type dt_in) : lap(lap_ref), dt(dt_in) {}
	template<typename grid_type, typename data_type, typename size_type>
	void operator()(size_type x, size_type y, grid_type & g, data_type Kd = 1.0)
	{
		g(x,y) += dt*Kd*lap(x,y);
	}
private:
	grid & lap;
	const grid::data_type dt;
};

class calc_lap //: public func2d
{
public:
	calc_lap(grid & lap_ref) : lap(lap_ref) {}
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
	{
		lap(x,y) = ( g(x-1,y) + g(x+1,y) + g(x,y-1) + g(x,y+1) - 4.0*g(x,y) );
	}
private:
	grid & lap;
};

class check //: public func2d
{
public:
	check(std::ostream & ostr_in, std::ostream & ostr2_in, std::ostream & ostr3_in) : ostr(ostr_in), ostr2(ostr2_in),
											  ostr3(ostr3_in), extended(true) {}
	check(std::ostream & ostr_in) : ostr(ostr_in), ostr2(std::cout), ostr3(std::cout), extended(false) {}
									//here cout is just a dummy.
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
	{
		grid::data_type wl = g.height();
		//TODO: revise to use region information instead of grid.
		ostr << y << ' ' << ( g(x,y) / (std::sin( (Pi/wl)*(y/*+1.0*/) ) ) ) << '\n';
		if(extended)
		{
			ostr2 << y << ' ' << g(x,y) << '\n';
			ostr3 << y << ' ' << std::sin( (Pi/wl)*(y+1.0) ) << '\n';
		}
	}
private:
	std::ostream & ostr, & ostr2, & ostr3;
	bool extended;
};

class freq //: public func2d
{
public:
	freq(std::ostream & ostr_in, grid::size_type n_points_in) : ostr(ostr_in), n_points(n_points_in), avg(0)
	{
		ostr << std::scientific << std::setprecision(grid::precision);
	}
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)	//prgrid::size_types out  "x  y  w"; updates average
	{
		grid::data_type w = -( (g(x-1,y) + g(x+1,y) + g(x,y-1) + g(x,y+1) - 4*g(x,y)) / g(x,y) );
		//w at this point is sometimes negative! (TE mode 0)
		if(w >= 0) //to prevent to printing of "nan" (because sqrt(-1) trumps even the checks within)
		{
            w = std::sqrt(w);
            ostr << x << ' ' << y << ' ' << w << '\n';
            if(is_number(w))
                avg += w / static_cast<grid::data_type>(n_points);
		}
	}
	template<typename T>
	bool is_number(const T & w)
	{ 
		return (w != std::numeric_limits<T>::quiet_NaN()) 
			&& (w == w) 
			&& w != std::numeric_limits<T>::signaling_NaN();
	}
	grid::data_type const average() { return avg; }

private:
	std::ostream & ostr;
	grid::size_type n_points;
	grid::data_type avg;
};

class to_stream //: public func2d
{
public:
	to_stream(std::ostream & ostr_in) : ostr(ostr_in)
	{
		ostr_in << std::scientific << std::setprecision(grid::precision);
	}
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
	{
	    ostr << x << '\t' << y << '\t' << g(x,y) << '\n';
	}
private:
	std::ostream & ostr;
};

class to_binary //: public func2d
{
public:
    to_binary(const std::string & path)
        : fout(path.c_str(), std::ios::binary),
          sizeofdatatype(static_cast<unsigned int>(sizeof(grid::data_type))),
          sizeofcoordcomponent(static_cast<unsigned short>(sizeof(unsigned short))),
          npertuple(2),
          badstate(false)
    {
        //The following 'if' condition may be problematic. Look closer before using binary save mode.
        if(!fout)
        {
            badstate = true;
            throw std::string("SHIT SHIT BAD PATH IN BINARY OUTPUT CONSTRUCTOR!\nTried ")+path;
        }
        else
        {
            fout.write( (char*)(&sizeofdatatype), sizeof(unsigned int) );
            fout.write( (char*)(&sizeofcoordcomponent), sizeof(unsigned short) );
            fout.write( (char*)(&npertuple), sizeof(unsigned short) );
        }
    }
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
    {
        if(badstate)
            throw std::string("bad func2d::to_binary attempted operator()");
        else
        {
            fout.write( (char*)(&x), sizeofcoordcomponent );
            fout.write( (char*)(&y), sizeofcoordcomponent );
            fout.write( (char*)(&(g(x,y))), sizeofdatatype);
        }
    }
    bool bad() const { return badstate; }
private:
    std::ofstream fout;
    const unsigned int sizeofdatatype;
    const unsigned short sizeofcoordcomponent;
    const unsigned short npertuple;
    bool badstate;
};

#ifdef CAVITY_USE_BOOST_RANDOM
class make_rand //: public func2d
{
public:
	make_rand(grid::size_type mx_in) : mx(mx_in), six(0,32767), die(rng, six)
	{
   		   //random integer [0 , 2^15-1]
	}
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
	{
		g(x,y) = (die() / 32767.0) * mx;
	}
private:
	grid::size_type mx;
	boost::mt19937 rng;
	boost::uniform_int<> six;
	boost::variate_generator<boost::mt19937&, boost::uniform_int<> > die;
};
#else
class make_rand //: public func2d
{
public:
    make_rand(grid::size_type mx_in) : mx(mx_in)
    {
    }
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
    {
        g(x,y) = static_cast<grid::data_type>(std::rand()) / RAND_MAX * mx;
        //g(x,y) -= mx/2;
    }
    grid::size_type mx;
};
#endif

class make_paraboloid //: public func2d
{
public:
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
	{
		g(x,y) = (x*y*(g.width() - x)*(g.height() - y)) / (g.height()*g.width());
	}
};

class make_dome //: public func2d
{
public:
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
	{
		g(x,y) = std::sqrt( std::abs(std::pow(g.width()/2.0,2.0) - std::pow((x - g.width()/2.0),2.0) - std::pow((y - g.height()/2.0),2.0)) );
	}
};

template<class unary_func, bool Kd_weighted = false>
class sum //: public func2d
{
public:
    sum(const unary_func & f_) : s(0.0), f(f_) {}
    sum() : s(0.0) {}
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
    {
        s += Kd_weighted ? Kd*f( g(x,y) ) : f( g(x,y) );
		//Perhaps use template specialization instead of this fork?
    }
    void reset() { s = 0.0; }
    const grid::data_type value() const { return s; }
private:
    grid::data_type s;
    unary_func f;
};

class inner_product //: public func2d
{
public:
    inner_product(grid & other) : g2(other), s(0) {}
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
    {
        s += g(x,y) * g2(x,y);
    }
    void reset()
    {
        s = 0.0;
    }
    const grid::data_type value() const { return s; }
private:
    grid & g2;
    grid::data_type s;
};

struct sub_mode //: public func2d
{
    sub_mode(grid & oldmode, grid::data_type sumofsq, grid::data_type inner_product)
        : old(oldmode), sosq(sumofsq), ip(inner_product) {}

    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
    {
    	g(x,y) -= (ip / sosq) * old(x,y);
    }
  private:
    grid & old;
    grid::data_type sosq;
	grid::data_type ip;
};

class max //: public func2d
{
public:
    max() : m(std::numeric_limits<grid::data_type>::min()) {}
    void reset() { m = std::numeric_limits<grid::data_type>::min(); }
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
    {
        if( g(x,y) > m ) m = g(x,y);
    }
    const grid::data_type & value() const { return m; }
private:
    grid::data_type m;
};

class min //: public func2d
{
public:
    min() : m(std::numeric_limits<grid::data_type>::max()) {}
    void reset() { m = std::numeric_limits<grid::data_type>::max(); }
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0)
    {
        if( g(x,y) < m ) m = g(x,y);
    }
    const grid::data_type & value() const { return m; }
private:
    grid::data_type m;
};

class set_to //: public func2d
{
public:
	set_to(grid::data_type val_in) : val(val_in) {}
	void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0) { g(x,y) = val; }
private:
	grid::data_type val;
};

class multiply //: public func2d
{
public:
    multiply(grid::data_type val_in) : val(val_in) {}
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0) { g(x,y) *= val; }
private:
	grid::data_type val;
};

class add //: public func2d
{
public:
    add(grid::data_type val_in) : val(val_in) {}
    void operator()(grid::size_type x, grid::size_type y, grid::base_type & g, grid::data_type Kd = 1.0) { g(x,y) += val; }
private:
	grid::data_type val;
};

} //namespace functions

#endif //GRID_H

